Clustered
bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset <- boostrap_state_time_group(data_subset,
sensitivity_anlys_post_tx_model_log_smoothed_time_subset,
coef_of_interest = rownames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,]),
nSim = 5000)
colnames(bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset) <-
rownames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,])
mean_coef_outcome_lin_post_tx_subset <- apply(bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset, 2, mean)
ci_coef_outcome_lin_post_tx_subset <- apply(bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset, 2, quantile, probs = c(.025, .975))
# write.csv(bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset,
# "./Data/all_clustered_bootstrap_smoothed_time_lin_tx_subset_3_11_22_nSim_5000.csv",
# row.names = FALSE)
ci_coef_outcome_lin_post_tx_subset <- cbind(t(ci_coef_outcome_lin_post_tx_subset),
mean_coef_outcome_lin_post_tx_subset)
colnames(ci_coef_outcome_lin_post_tx_subset) <- c("conf.low", "conf.high", "estimate")
ci_coef_outcome_lin_post_tx_subset <- data.frame(ci_coef_outcome_lin_post_tx_subset)
ci_coef_outcome_lin_post_tx_subset$term <- rownames(ci_coef_outcome_lin_post_tx_subset)
ci_coef_outcome_lin_post_tx_subset$ci_95 <- paste("95% CI = (",
format(round(ci_coef_outcome_lin_post_tx_subset$conf.low, 3), nsmall = 3),
",",
format(round(ci_coef_outcome_lin_post_tx_subset$conf.high,3), nsmall = 3),
")", sep = "")
# write.csv(ci_coef_outcome_lin_post_tx_subset,
# "./Data/clustered_bootstrap_smoothed_time_lin_tx_subset_3_11_22_nSim_5000.csv",
# row.names = FALSE)
bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset <-
read.csv("./Data/all_clustered_bootstrap_smoothed_time_lin_tx_subset_3_11_22_nSim_5000.csv")
#first examine the distribution of coefficients
bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset_long <- bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset %>%
pivot_longer(cols = everything(),
names_to = "policy",
values_to = "coef_est")
ggplot(bootstrap_smoothed_eff_log_outcome_lin_post_tx_subset_long, aes(y = coef_est, x = policy)) +
geom_boxplot() +
labs(x = "Policies", y = "Coefficient Estimates",
title = "Clustered Bootstrap Distribution of Coefficient Estimates",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
axis.text.x = element_text(size = 5, angle = 45))

coef_w_ci_smoothed_time_log_outcome_lin_post_tx_subset <-
read.csv("./Data/clustered_bootstrap_smoothed_time_lin_tx_subset_3_11_22_nSim_5000.csv",
row.names = "term")
coef_w_ci_smoothed_time_log_outcome_lin_post_tx_subset$term <- rownames(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_subset)
#multiply coefficients by 10 if they're the linear effects
ten_pd_eff_subset <- coef_w_ci_smoothed_time_log_outcome_lin_post_tx_subset[
rownames(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_subset) %in%
c("lag_num_pd_w_naloxone_yes",
"lag_num_pd_w_naloxone_no",
"lag_num_pd_w_med_marijuana",
"lag_num_pd_w_rec_marijuana",
"lag_num_pd_w_gsl",
"lag_num_pd_w_pdmp",
"lag_num_pd_w_medicaid",
"lag_num_pd_w_tx"),]
ten_pd_eff_subset$estimate <- ten_pd_eff_subset$estimate*10
ten_pd_eff_subset$conf.high <- ten_pd_eff_subset$conf.high*10
ten_pd_eff_subset$conf.low <- ten_pd_eff_subset$conf.low*10
ten_pd_eff_subset$ci_95 <- paste("95% CI = (", format(round(ten_pd_eff_subset$conf.low, 3), nsmall = 3), ", ",
format(round(ten_pd_eff_subset$conf.high, 3), nsmall = 3), ")", sep = "")
rownames(ten_pd_eff_subset) <- ten_pd_eff_subset$term <- c("num_pd_w_naloxone_yes_per_5_yr",
"num_pd_w_naloxone_no_per_5_yr",
"num_pd_w_med_marijuana_per_5_yr",
"num_pd_w_rec_marijuana_per_5_yr",
"num_pd_w_gsl_per_5_yr",
"num_pd_w_pdmp_per_5_yr",
"num_pd_w_medicaid_per_5_yr",
"num_pd_w_tx_per_5_yr")
combined_ten_pd_df_subset <- rbind(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_subset, ten_pd_eff_subset)
combined_ten_pd_df_subset_to_combine <- combined_ten_pd_df_subset %>%
filter(!(term %in% c("lag_num_pd_w_naloxone_yes",
"lag_num_pd_w_naloxone_no",
"lag_num_pd_w_med_marijuana",
"lag_num_pd_w_rec_marijuana",
"lag_num_pd_w_gsl",
"lag_num_pd_w_pdmp",
"lag_num_pd_w_medicaid",
"lag_num_pd_w_tx")))
#put model estimates in the plot
model_effect_estimates_subset <- sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,]
model_effect_estimates_subset_to_combine <- model_effect_estimates_subset[
!rownames(model_effect_estimates_subset) %in%
c("lag_num_pd_w_naloxone_yes",
"lag_num_pd_w_naloxone_no",
"lag_num_pd_w_med_marijuana",
"lag_num_pd_w_rec_marijuana",
"lag_num_pd_w_gsl",
"lag_num_pd_w_pdmp",
"lag_num_pd_w_medicaid",
"lag_num_pd_w_tx"),]
combined_ten_pd_model_est_subset <- sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[
rownames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset) %in%
c("lag_num_pd_w_naloxone_yes",
"lag_num_pd_w_naloxone_no",
"lag_num_pd_w_med_marijuana",
"lag_num_pd_w_rec_marijuana",
"lag_num_pd_w_gsl",
"lag_num_pd_w_pdmp",
"lag_num_pd_w_medicaid",
"lag_num_pd_w_tx"),]
rownames(combined_ten_pd_model_est_subset) <- combined_ten_pd_model_est_subset$term <- c("num_pd_w_naloxone_yes_per_5_yr",
"num_pd_w_naloxone_no_per_5_yr",
"num_pd_w_med_marijuana_per_5_yr",
"num_pd_w_rec_marijuana_per_5_yr",
"num_pd_w_gsl_per_5_yr",
"num_pd_w_pdmp_per_5_yr",
"num_pd_w_medicaid_per_5_yr",
"num_pd_w_tx_per_5_yr")
combined_ten_pd_model_est_subset$estimate <- combined_ten_pd_model_est_subset$estimate*10
combined_ten_pd_model_est_subset$conf.low <- combined_ten_pd_model_est_subset$conf.low*10
combined_ten_pd_model_est_subset$conf.high <- combined_ten_pd_model_est_subset$conf.high*10
combined_ten_pd_model_est_subset$ci_95 <- paste("95% CI = (",
format(round(combined_ten_pd_model_est_subset$conf.low, 3), nsmall = 3),
", ",
format(round(combined_ten_pd_model_est_subset$conf.high, 3), nsmall = 3), ")", sep = "")
model_effect_estimates_subset_to_combine <- rbind(model_effect_estimates_subset_to_combine, combined_ten_pd_model_est_subset)
ggplot() +
geom_point(combined_ten_pd_df_subset_to_combine, mapping = aes(x = estimate, y = 16:1, color = "Bootstrap Estimate")) +
geom_linerange(combined_ten_pd_df_subset_to_combine,
mapping = aes(xmin = conf.low, xmax = conf.high, y = 16:1, color = "Bootstrap Estimate")) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(color = "black"),
legend.position = "bottom") +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
geom_point(model_effect_estimates_subset_to_combine,
mapping = aes(x = estimate, y = 16:1, color = "Model Estimate")) +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Model Coefficient and Bootstrap 95% CI of
Model Using Smoothed Time Effects,
Linear Policy Effects,
Using Clustered Bootstrap",
color = "Estimate Type") +
# scale_color_grey() +
geom_text(model_effect_estimates_subset_to_combine,
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 1.5, y = 16:1), size = 3) +
geom_text(combined_ten_pd_df_subset_to_combine,
mapping = aes(label = ci_95, x = 2.6, y = 16:1), size = 3) +
xlim(-2, 3.3) +
scale_y_continuous(breaks = 1:16, labels = combined_ten_pd_df_subset_to_combine$term[16:1])

ggplot() +
# geom_point(combined_ten_pd_df_subset, mapping = aes(x = estimate, y = 16:1, color = "Bootstrap Estimate")) +
geom_linerange(combined_ten_pd_df_subset_to_combine, mapping = aes(xmin = conf.low, xmax = conf.high, y = 16:1)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(color = "black"),
legend.position = "bottom") +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
geom_point(model_effect_estimates_subset_to_combine,
mapping = aes(x = estimate, y = 16:1)) +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Model Coefficient and Bootstrap 95% CI of
Model Using Smoothed Time Effects,
Linear Policy Effects,
Using Clustered Bootstrap for Subset Data",
color = "Estimate Type") +
# scale_color_grey() +
geom_text(model_effect_estimates_subset_to_combine,
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 1.5, y = 16:1), size = 3) +
geom_text(combined_ten_pd_df_subset_to_combine,
mapping = aes(label = ci_95, x = 2.6, y = 16:1), size = 3) +
xlim(-2, 3.3) +
scale_y_continuous(breaks = 1:16, labels = combined_ten_pd_df_subset_to_combine$term[16:1])

ggplot() +
# geom_point(combined_ten_pd_df_subset, mapping = aes(x = estimate, y = 16:1, color = "Bootstrap Estimate")) +
geom_linerange(combined_ten_pd_df_subset_to_combine,
mapping = aes(xmin = conf.low, xmax = conf.high, y = 16:1, color = "Subset Data"),
alpha = 0.5) +
geom_linerange(model_effect_estimates_to_be_combined,
mapping = aes(xmin = conf.low, xmax = conf.high, y = 16:1, color = "Full Data"), alpha = 0.5) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(color = "black"),
legend.position = "bottom") +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
geom_point(model_effect_estimates_subset_to_combine,
mapping = aes(x = estimate, y = 16:1, color = "Subset Data")) +
geom_point(model_effect_estimates_to_be_combined,
mapping = aes(x = estimate, y = 16:1, color = "Full Data")) +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Model Coefficient and Bootstrap 95% CI of
Model Using Smoothed Time Effects,
Linear Policy Effects,
Using Clustered Bootstrap for Full and Subset Data",
color = "Estimate Type") +
# scale_color_grey() +
# geom_text(model_effect_estimates_subset,
# mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 1.5, y = 16:1), size = 3) +
# geom_text(combined_ten_pd_df_subset,
# mapping = aes(label = ci_95, x = 2.6, y = 16:1), size = 3) +
# xlim(-2, 3.2) +
scale_y_continuous(breaks = 1:16, labels = combined_ten_pd_df_subset_to_combine$term[16:1])

sd_for_attr_death_subset <- merge(model_effect_estimates_subset[,c("term", "estimate")],
combined_ten_pd_df_subset[,c("conf.low", "term", "conf.high", "ci_95")],
by = "term")
rownames(sd_for_attr_death_subset) <- sd_for_attr_death_subset$term
bootstrap_attributable_deaths_subset <- attr_death_compute(data_subset,
sd_for_attr_death_subset,
lin_model = TRUE,
tx_name = c("Intervention_Redefined", "lag_num_pd_w_tx"))
bootstrap_attributable_deaths_subset <- merge(bootstrap_attributable_deaths_subset, date_data_subset,
by.x = "Time_Period", by.y = "Time_Period_ID")
bootstrap_attributable_deaths_subset <- bootstrap_attributable_deaths_subset %>%
group_by(year = year(Time_Period_Start)) %>%
summarise(lives_saved = sum(attr_deaths),
lives_saved_lb = sum(attr_deaths_lb),
lives_saved_ub = sum(attr_deaths_ub))
ggplot(bootstrap_attributable_deaths_subset, aes(x = year)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = lives_saved, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = lives_saved_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = lives_saved_ub, linetype = "95% CI")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Per Year
Using Model with Smoothed Time Effects,
Log Probability of Drug Overdose Death,
Estimated Using Clustered Bootstrap on Subset Data",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))

ggplot() +
geom_line(bootstrap_attributable_deaths_subset,
mapping = aes(x = year, y = lives_saved, linetype = "Estimate", color = "Subset Data")) +
geom_line(bootstrap_attributable_deaths_subset,
mapping = aes(x = year, y = lives_saved_lb, linetype = "95% CI", color = "Subset Data")) +
geom_line(bootstrap_attributable_deaths_subset,
mapping = aes(x = year, y = lives_saved_ub, linetype = "95% CI", color = "Subset Data")) +
geom_line(bootstrap_attributable_deaths_summary[bootstrap_attributable_deaths_summary$year <= 2014, ],
mapping = aes(x = year, y = total_lives_saved, linetype = "Estimate", color = "Full Data")) +
geom_line(bootstrap_attributable_deaths_summary[bootstrap_attributable_deaths_summary$year <= 2014, ],
mapping = aes(x = year, y = total_lives_saved_lb, linetype = "95% CI", color = "Full Data")) +
geom_line(bootstrap_attributable_deaths_summary[bootstrap_attributable_deaths_summary$year <= 2014, ],
mapping = aes(x = year, y = total_lives_saved_ub, linetype = "95% CI", color = "Full Data")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Per Year
Using Model with Smoothed Time Effects
and Linear Policy Effects,
Estimated Using Clustered Bootstrap
for Subset Data and Full Data until 2015",
linetype = "", color = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))

Unit Level
bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset <- boostrap_state_time_unit(data_subset,
sensitivity_anlys_post_tx_model_log_smoothed_time_subset,
coef_of_interest = rownames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,]))
colnames(bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset) <-
rownames(sensitivity_anlys_post_tx_sd_and_ci_log_smoothed_time_subset[51:66,])
mean_coef_outcome_lin_post_tx_unit_subset <- apply(bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset, 2, mean)
ci_coef_outcome_lin_post_tx_unit_subset <- apply(bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset, 2,
quantile, probs = c(.025, .975))
summary_coef_outcome_lin_post_tx_unit_subset <- apply(bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset, 2,
summary)
# write.csv(bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset ,
# "./Data/all_unit_bootstrap_smoothed_time_lin_tx_subset_3_10_22.csv",
# row.names = FALSE)
coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset <- cbind(t(ci_coef_outcome_lin_post_tx_unit_subset),
mean_coef_outcome_lin_post_tx_unit_subset)
colnames(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset) <- c("conf.low", "conf.high", "estimate")
coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset <- data.frame(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset)
coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset$term <- rownames(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset)
coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset$ci_95 <- paste("95% CI = (",
format(round(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset$conf.low,
3), nsmall = 3),
",",
format(round(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset$conf.high,3), nsmall = 3),
")", sep = "")
# write.csv(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset ,
# "./Data/unit_bootstrap_smoothed_time_lin_tx_subset_3_10_22.csv",
# row.names = FALSE)
bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset <- read.csv("./Data/all_unit_bootstrap_smoothed_time_lin_tx_subset_3_10_22.csv")
#first examine the disribution of coefficients
bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset_long <- bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset %>%
pivot_longer(cols = everything(),
names_to = "policy",
values_to = "coef_est")
ggplot(bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset_long, aes(y = coef_est, x = policy)) +
geom_boxplot() +
labs(x = "Policies", y = "Coefficient Estimates",
title = "Bootstrap Distribution of Coefficient Estimates",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
axis.text.x = element_text(size = 5, angle = 45))

ggplot(bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset_long[
!bootstrap_smoothed_eff_log_outcome_lin_post_tx_unit_subset_long$policy %in%
c("Recreational_Marijuana_Redefined", "lag_num_pd_w_rec_marijuana"),], aes(y = coef_est, x = policy)) +
geom_boxplot() +
labs(x = "Policies", y = "Coefficient Estimates",
title = "Bootstrap Distribution of Coefficient Estimates Without Recreational Marijuana Policies",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
axis.text.x = element_text(size = 5, angle = 45))

coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset <- read.csv("./Data/unit_bootstrap_smoothed_time_lin_tx_subset_3_10_22.csv",
row.names = "term")
coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset$term <- rownames(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset)
dwplot(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset, colour = "black") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_vline(aes(xintercept = 0), linetype = "dashed") +
labs(y = "Term", x = "Coefficients and 95% Confidence Intervals",
title = "Coefficient and 95% CI of Linear Post Tx Effects,
Smoothed Time Effects,
Using Unit Bootstrap, Subset Data") +
scale_color_grey() +
geom_text(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset,
mapping = aes(label = format(round(estimate, 3), nsmall = 3), x = 0.55, y = 16:1), size = 3) +
geom_text(coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset,
mapping = aes(label = ci_95, x = 1, y = 16:1), size = 3) +
xlim(-.65, 1.2)

bootstrap_attributable_deaths_unit_subset <- attr_death_compute(data_subset,
coef_w_ci_smoothed_time_log_outcome_lin_post_tx_unit_subset,
lin_model = TRUE,
tx_name = c("Intervention_Redefined", "lag_num_pd_w_tx"))
bootstrap_attributable_deaths_unit_subset <- merge(bootstrap_attributable_deaths_unit_subset, date_data_subset,
by.x = "Time_Period", by.y = "Time_Period_ID")
ggplot(bootstrap_attributable_deaths_unit_subset, aes(x = Time_Period_Start)) +
# geom_point(aes(y = attr_deaths)) +
geom_line(aes(y = attr_deaths, linetype = "Estimate")) +
# geom_point(aes(y = attr_deaths_lb)) +
geom_line(aes(y = attr_deaths_lb, linetype = "95% CI")) +
# geom_point(aes(y = attr_deaths_ub)) +
geom_line(aes(y = attr_deaths_ub, linetype = "95% CI")) +
labs(x = "Date", y = "Lives Saved",
title = "Estimated Number of Lives Saved Using Smoothed Time Effects,
Log Probability of Drug Overdose Death,
Estimated Using Unit Bootstrap",
linetype = "") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_linetype_manual(values = c("dashed", "solid"))
